program test
  use mod_bivar
  use mod_workspace
  USE DFPORT 
  implicit none

integer,parameter::tndp=100*90
  real(8),dimension(tndp):: xd,yd,zd,zd2
  integer::i,j

  ndp =  tndp
  nxi = 100
  nyi = 130

  xd(1)=rand(789)
  do i = 1, ndp
     xd(i)=4*rand()-2
     yd(i)=4*rand()-2
  enddo
  zd =  xd*exp(-xd**2-yd**2)
!write(9,*)xd,yd,zd
  call bivar_init1()

  do i = 1, nxi
     xi(i) = 4*real ( i - 1, kind = 8 ) / real ( nxi - 1, kind = 8 )-2
  end do

  do j = 1, nyi
     yi(j) = 4*real ( j - 1, kind = 8 ) / real ( nyi - 1, kind = 8 )-2
  end do
  write(*,*)"griddata1:"
  call timestamp()
  call griddata1(xd,yd,zd,zi)
  call timestamp()

  call timestamp()
  call griddata1(xd,yd,zd,zi)
  call timestamp()

 ! call bivar_init2()
 ! write(*,*)minval(xi),minval(xd),maxval(xi),maxval(xd)
 ! write(*,*)minval(yi),minval(yd),maxval(yi),maxval(yd)

 ! write(*,*)"griddata2:"
 ! call timestamp()
 ! call griddata2(xd,yd,zd2,zi)
 ! call timestamp()
 ! zd2=zd-zd2

 ! write(*,*)"error= " ,maxval(zd2)
 ! write(*,*)"griddata2, second time:"
 ! call timestamp()
 ! call griddata2(xd,yd,zd2,zi)
 ! call timestamp()
  



  write(*,*)"bicubic:"
  call timestamp()
  call init_interp2(nxi,nyi,xi,yi,tndp,xd,yd)
  do i=1,100
  call interp2(zi,zd2)
  enddo
  call timestamp()
  !write(10,*)zd2
  zd2=zd-zd2
  write(*,*)"error= " ,minval(zd2)
  

  !call delete_bivar_workspace()
    call timestamp()
    do i=1,100
  call imslinterp2(nxi,xi,nyi,yi,zi,ndp,xd,yd,zd2)
    enddo
	  call timestamp()
  zd2=zd-zd2
  write(*,*)"error= " ,minval(zd2)



end program test



!------------------------------------------------
subroutine imslinterp2(nxdata,xdata,nydata,ydata,fvalue,ndp,xd,yd,zd)
!------------------------------------------------
use imsl
implicit none
integer,parameter :: kxord=5,kyord=5,maxsize=1000
integer:: nxdata,nydata,nxvec,ndp,i
real*8,dimension(maxsize)::xknot,yknot
real*8,dimension(maxsize,maxsize)::bscoef
real*8::xdata(nxdata),ydata(nydata),fvalue(nxdata,nydata),xd(ndp),yd(ndp),zd(ndp)


!CALL ERSET (0 , 1, 0)
call dbsnak(nxdata, xdata, kxord, xknot) 
call dbsnak(nydata, ydata, kyord, yknot)
call dbs2in (nxdata, xdata, nydata, ydata, fvalue, nxdata, kxord,&
kyord, xknot, yknot,bscoef)

!ierr=IERCD()
!if (ierr .gt. 0) then
! pause "error"
!endif
do i=1,ndp
    zd(i)= dbs2vl (xd(i), yd(i), kxord, kyord, xknot,&
    yknot, nxdata, nydata, bscoef)  
enddo


end subroutine imslinterp2

